#include "precompiled.h"
////
//  Standard headers
// 
#include <iterator>   // for stream_iterators, inserter
#include <iostream>
	using std::cout;
	using std::cerr;
	using std::cin;
	using std::ofstream;
	using std::ifstream;
#include <string>
	using std::string;
#include <vector>
	using std::vector;
#include <map>
	using std::map;
#include <utility>
#include <algorithm>
	using std::min;
	using std::max;
#include <numeric>
#include <stdexcept>
	using std::runtime_error;
#include <boost/range/algorithm_ext/push_back.hpp>
    using boost::push_back;


//
//  My common headers
// 
#include <progress_indicator.h>
#include <commify.h>
#include <glob_files.h>
#include <open_stream_name_with_substitutions.h>


//
//  Boost headers
// 
#include <boost/filesystem/path.hpp>
	using boost::filesystem::path;
#include <boost/filesystem/convenience.hpp>
	using boost::filesystem::basename;
	using boost::filesystem::is_directory;
	using boost::filesystem::create_directories;
#include <boost/regex.hpp>
	using boost::regex;
	using boost::regex_search;
	using boost::match_results;

#include <boost/algorithm/string/join.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/foreach.hpp>
#include <boost/tr1/tuple.hpp>
#include <boost/tr1/unordered_map.hpp>
	using std::tr1::unordered_map;
#include <boost/unordered_set.hpp>
	using boost::unordered_set;
#include <boost/range/algorithm.hpp>
	using boost::sort;
#include <boost/format.hpp>
		using boost::format;
		using boost::io::group;
#include <boost/algorithm/string/replace.hpp>
		using boost::algorithm::replace_all;
#include <boost/iostreams/device/file.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/mapped_file.hpp>
		namespace io = boost::iostreams;

#include <boost/assign/list_of.hpp> // for 'list_of()'
		using namespace boost::assign; // bring 'list_of()' into scope

#include <boost/assert.hpp> 





//
//  tr1 headers
// 
#include <boost/tr1/tuple.hpp>
	using std::tr1::tuple;
	using std::tr1::make_tuple;
	using std::tr1::get;
	

#include "intervals_per_contig.h"
	using intervals_per_contig::t_interval;
	using intervals_per_contig::t_intervals;
	using intervals_per_contig::t_intervals_per_contig;
	using intervals_per_contig::t_vec_intervals;
	using intervals_per_contig::t_map_str_intervals;

#include "read_options.h"




//!_______________________________________________________________________________________
//!     
//!		read_intervals_from_file
//!     
//!		\author lg (1/26/2011)
//!     
//!     \param opt 
//!     \param input_file_name 
//!     \param intervals_per_contig 
//!     \param zero_based_input 
//!_______________________________________________________________________________________
void read_intervals_from_file(	prog_options&			opt, 
								const string&			input_file_name, 
								t_intervals_per_contig& intervals_per_contig,
								bool					zero_based_input)
{
	io::stream<io::mapped_file_source> f_input(input_file_name);
	if (!f_input)
		throw std::logic_error("Can't open " + input_file_name);
	read_intervals(	f_input, intervals_per_contig, zero_based_input);
	intervals_per_contig.combine_overlapping ();
}




//!_______________________________________________________________________________________
//!     
//!		get_contig_sizes
//! 
//! 		returns empty map if no file specified
//!     
//!		\author lg (1/24/2011)
//!     
//!     \param opt 
//!     
//!     \return unordered_map<string,unsigned> 
//!_______________________________________________________________________________________
unordered_map<string, unsigned> get_contig_sizes(prog_options& opt)
{
	unordered_map<string, unsigned> contig_sizes;

	string file_name = opt.contig_sizes_file_name;
	if (!file_name.length())
	{
		opt.verbose(3) << "No contig file sizes specified.\n";
		return contig_sizes;
	}
	opt.both(3) << "Reading contig sizes  from " << file_name << "...\n";
	opt.both(3) << "  !!Assumes zero based coordinates!!\n";

	t_intervals_per_contig	intervals_per_contig;
	read_intervals_from_file(	opt, 
								file_name, 
								intervals_per_contig, true);

	vector<string> contig_names = intervals_per_contig.get_sorted_contigs();
	unsigned cnt_sizes = intervals_per_contig.size();
	if (cnt_sizes != contig_names.size())
		opt.both() << "WARNING! Some contigs have more than one size\n";

	//	only keep 1-99,X
	opt.both(3) << "    " << contig_names.size() << " contigs\n";
	opt.log() << "    contigs = [ " << boost::join(contig_names, ", ") << " ]\n";


	// 
	// only take the first size
	// 
	BOOST_FOREACH(const string& contig_name, contig_names) 
	{
		t_interval interval = intervals_per_contig.data[contig_name].data[0];
		contig_sizes[contig_name] = interval.end - interval.beg;
	}
	return contig_sizes;
}





//!_______________________________________________________________________________________
//!     
//!		summarise_intervals
//! 
//! 		returns summary (total) string as well as printing out to log
//! 
//!		\author lg (1/26/2011)
//!     
//!     \param opt 
//!     \param interval_set_name 
//!     \param intervals_per_contig 
//!     \param genome_contig_sizes 
//!     \param total_verbosity 
//!     
//!     \return string 
//!_______________________________________________________________________________________

string summarise_intervals(	prog_options&							opt, 
							const string&							interval_set_name,
							unsigned 								max_interval_set_name_len,
							t_intervals_per_contig&					intervals_per_contig, 
							const unordered_map<string, unsigned>&	genome_contig_sizes,
							unsigned								total_verbosity)
{
	opt.both(3) << "\n" << string (80, '_') << "\n" << interval_set_name << "\n";
	vector<string> contig_names = intervals_per_contig.get_sorted_contigs();
	opt.both(3)	<< contig_names.size() << " contigs\n";
	opt.log()	<< "    contigs = [ " << boost::join(contig_names, ", ") << " ]\n";

	//
	//	get longest contig name
	std::size_t longest_contig_name_len = 0;
	BOOST_FOREACH(const string& contig_name, contig_names) 
		longest_contig_name_len = std::max(longest_contig_name_len, contig_name.size());
	longest_contig_name_len = std::max(longest_contig_name_len, interval_set_name.length());
	longest_contig_name_len = std::max(longest_contig_name_len, std::size_t(max_interval_set_name_len));
	longest_contig_name_len = std::max(longest_contig_name_len, std::size_t(20));
	longest_contig_name_len +=1;

	string contig_name_format_str = (format("  %%1$-%1$ds:") % longest_contig_name_len).str();
	 
	// 
	// running total for contigs with intervals and contig size
	// 
	unsigned total_size_with_contigs_size = 0;
	unsigned total_genome_size_with_intervals  = 0;
	BOOST_FOREACH(const string& contig_name, contig_names) 
	{
		unsigned contig_size	= intervals_per_contig.data[contig_name].len();
		unsigned contig_cnt		= intervals_per_contig.data[contig_name].size();
		if (genome_contig_sizes.count(contig_name))
		{
			unsigned genome_contig_size = genome_contig_sizes.find(contig_name)->second;
			opt.both(3) << format(contig_name_format_str + " %2$13s intervals %3$13s / %4$13s bp = %5$5.2f%%\n")
								% contig_name 
								% commify(contig_cnt  						)
								% commify(contig_size 						)
								% commify(genome_contig_size				)
								% (contig_size * 100.0 / genome_contig_size	);
			total_size_with_contigs_size		+= contig_size;
			total_genome_size_with_intervals	+= genome_contig_size;

		}
		else
			opt.both(3) << format(contig_name_format_str + " %2$13s intervals %3$13s bp\n")
								% contig_name 
								% commify(contig_cnt  						)
								% commify(contig_size 						);
	}

	// 
	// totals for all contigs
	// 
	unsigned total_size	= intervals_per_contig.len();
	unsigned total_cnt	= intervals_per_contig.size();
	
	string total_string;
	if (total_genome_size_with_intervals)
	{
		if (total_size_with_contigs_size != total_size)
		{
			opt.both(3)  << format(contig_name_format_str + " %2$13s           %3$13s  %4$13s       %5$13s      %6$s\n")
								% string(5, '-') % string(7, '-') % string(13, '-') % string(13, '-') % string(13, '-') % string(6, '-');
			total_string = (format(contig_name_format_str + " %2$13s intervals %3$13s (%4$13s bp) / %5$13s bp = %6$5.2f%%\n")
								% interval_set_name
								% commify(total_cnt  						)
								% commify(total_size 						)
								% commify(total_size_with_contigs_size		)
								% commify(total_genome_size_with_intervals	)	
								% (total_size_with_contigs_size * 100.0 / total_genome_size_with_intervals)).str();
			opt.both(total_verbosity)  << total_string;
		}
		else
		{
			opt.both(3)  << format(contig_name_format_str + " %2$13s           %3$13s   %4$13s      %5$s\n")
								% string(5, '-') % string(7, '-') % string(13, '-') % string(13, '-') % string(6, '-');
			total_string = (format(contig_name_format_str + " %2$13s intervals %3$13s / %4$13s bp = %5$5.2f%%\n")
								% interval_set_name
								% commify(total_cnt  						)
								% commify(total_size 						)
								% commify(total_genome_size_with_intervals	)	
								% (total_size_with_contigs_size * 100.0 / total_genome_size_with_intervals)).str();
			opt.both(total_verbosity)  << total_string;
		}
	}
	else
	{
		opt.both(3)  << format(contig_name_format_str + " %2$13s           %3$13s\n")
							% string(5, '-') % string(7, '-') % string(13, '-');
		total_string = (format(contig_name_format_str + " %2$13s intervals %3$13s bp\n")
							% interval_set_name
							% commify(total_cnt  						)
							% commify(total_size 						)).str();
		opt.both(total_verbosity)  << total_string;
	}
	opt.both(3) << string (80, '_') << "\n";
	return total_string;
}


//!_______________________________________________________________________________________
//!     
//!		get_interval_set_names
//!			
//! 		uses either regular expression or filename for the interval set name
//!     
//!		\author lg (1/26/2011)
//!     
//!     \param opt 
//!     \param input_file_name 
//!     \param interval_set_names 
//!     
//!     \return unsigned max interval set name length
//!_______________________________________________________________________________________

unsigned get_interval_set_names(prog_options&						opt, 
								const vector<string>& 				input_file_names,
								vector<string>& 					interval_set_names)
{
	// 
	//	get interval set name 
	// 
	std::size_t max_interval_set_name_len = 0U;
	if (opt.name_regex.length())
	{
		regex name_regex(opt.name_regex);
		BOOST_FOREACH(const string& input_file_name, input_file_names) 
		{
			string interval_set_name = boost::regex_replace(input_file_name, name_regex, 
															opt.name_sub.length() ? opt.name_sub : "$&",
															boost::match_default | boost::format_no_copy | 
															boost::match_perl);
			if (interval_set_name.length() == 0)
				throw std::logic_error((format("The regular expression %1$s does not match the input file name [%2$s]")
										% (opt.name_sub.length() ? 
											("'" + opt.name_regex + "' / '" + opt.name_sub + "'") :  
											("'" + opt.name_regex + "'")								)
										% input_file_name).str());
			max_interval_set_name_len = max(max_interval_set_name_len, interval_set_name.size());
			interval_set_names.push_back(interval_set_name);
		}
	}
	else if (opt.name_sub.length())
	{
		max_interval_set_name_len = opt.name_sub.length();
		if (input_file_names.size() > 1)
			opt.both() << "WARNING! All interval sets have the same name.\n";
		for (unsigned i = 0; i < input_file_names.size(); ++i)
			interval_set_names.push_back(opt.name_sub);
	}
	else
		BOOST_FOREACH(const string& input_file_name, input_file_names) 
		{
			string interval_set_name = basename(input_file_name);	
			max_interval_set_name_len = max(max_interval_set_name_len, interval_set_name.size());
			interval_set_names.push_back(interval_set_name);
		}
	return max_interval_set_name_len;
}

void read_second_intervals_per_contig(	prog_options&							opt, 
										const string&							input_file_name,
										unsigned 								max_interval_set_name_len,
										const unordered_map<string, unsigned>&	genome_contig_sizes,
										string& 								input_file_name2,
										t_intervals_per_contig&					intervals_per_contig2,
										const map<string, string>& 				custom						)
{
	//
	//  only read in intervals if different file name
	// 
	string new_input_file_name2 = get_file_name_substitutions(	opt.input_file_name2,
																input_file_name, custom);
	if (input_file_name2 == new_input_file_name2)
		return;
	input_file_name2 = new_input_file_name2;

	// read intervals
	intervals_per_contig2.clear();
	read_intervals_from_file(opt, input_file_name2, intervals_per_contig2, opt.zero_based_input2);

	// 
	// summarise new file
	// 
	string interval_set_name2 = basename(input_file_name2);
	summarise_intervals(opt, interval_set_name2, max_interval_set_name_len, intervals_per_contig2, genome_contig_sizes, 2);

}

string get_operation_string(prog_options& opt)
{
	switch (opt.operation)
	{
	case ePASS_THROUGH:
	case eSUMMARISE:
		return string();
	case eEXPAND:
		return (format(" (Expanded by %1$d)") % opt.offset).str();
	case eBRACKET:
		return (format(" (Bracketed by %1$d)") % opt.offset).str();
	case eUNION:
		return " (Union)";
	case eINTERSECT:
		return " (Intersection)";
	case eSUBTRACT:
		return " (Diff)";
	case eCLOSE_NEIGHBOURS:
		return (format(" (Join < %1$d)") % opt.offset).str();
	case eREMOVE_SHORT:
		return (format(" (Remove < %1$d)") % opt.offset).str();
	}
	throw std::logic_error("Unknown operation.");
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// 
//      Main
// 
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888 
const unsigned cnt_lines_per_chunk = 1500;
int main (int argc, char *argv[])
{
	try
	{
		// 
		//  Read options
		//  	return if user chose help or version
		// 
		prog_options opt;
		if (!read_options(opt, argc, argv))
			return 0;

		string input_file_name2;
		t_intervals_per_contig intervals_per_contig2;

		//
		//  get file names for input intervals 
		//
		vector<string> input_file_names = glob_files::glob(opt.input_file_name);
		if (input_file_names.size() == 0)
			throw std::logic_error("No files matching " + opt.input_file_name);

		//
		//	Get interval set names
		// 
		vector<string> interval_set_names;
		unsigned max_interval_set_name_len = get_interval_set_names(opt, input_file_names, interval_set_names);

		string operation_str = get_operation_string(opt);
		max_interval_set_name_len += operation_str.length();


		string summary_string = "\n\n" + string (80, '*') + "\nSummary:\n";
		for (unsigned ii = 0; ii < input_file_names.size(); ++ii)
		{
			const string& input_file_name 	= input_file_names[ii];
			const string& interval_set_name	= interval_set_names[ii];
			// 
			//	read intervals 
			// 
			t_intervals_per_contig intervals_per_contig ;
			opt.both(3) << "Reading intervals from " << input_file_name << "...\n";
			read_intervals_from_file(opt, input_file_name, intervals_per_contig, opt.zero_based_input);

			//	
			// Get contig sizes 
			// 
			unordered_map<string, unsigned> genome_contig_sizes = get_contig_sizes(opt);


			//
			//	Prepare for file name substitution for output and second input file name
			// 
			map<string, string> custom;
			custom["[NAME]"] = interval_set_name;

			// 
			//	input_file_name2 
			// 
			if (opt.operation >= eUNION)
				read_second_intervals_per_contig(opt, input_file_name, max_interval_set_name_len,
												genome_contig_sizes, 
												input_file_name2, intervals_per_contig2, custom);
			
			//
			//	summarise intervals before operation 
			// 
			if (opt.operation != eSUMMARISE && opt.operation != ePASS_THROUGH)
				summarise_intervals(opt, basename(input_file_name), max_interval_set_name_len, intervals_per_contig, genome_contig_sizes, 2);


			// do operation
			switch (opt.operation)
			{
			case ePASS_THROUGH:
			case eSUMMARISE:
				break;
			//case eCONTRACT:
			//	break;
			case eEXPAND:
				intervals_per_contig.expand_by(opt.offset);
				break;
			case eBRACKET:
				intervals_per_contig.bracket_by(opt.offset);
				break;
			case eUNION:
				intervals_per_contig.combine_with(intervals_per_contig2);
				break;
			case eINTERSECT:
				intervals_per_contig.intersect_with(intervals_per_contig2);
				break;
			case eSUBTRACT:
				intervals_per_contig.subtract(intervals_per_contig2);
				break;
			case eCLOSE_NEIGHBOURS:
				intervals_per_contig.combine_close_neighbours(opt.offset);
				break;
			case eREMOVE_SHORT:
				intervals_per_contig.remove_intervals_smaller_than(opt.offset);
				break;
			}



			//
			//	summarise intervals after operation 
			// 
			string transformed_interval_set_name = interval_set_name + operation_str;
			unsigned spacer_len = max_interval_set_name_len - (interval_set_name + operation_str).length();
			if (spacer_len)
				transformed_interval_set_name = interval_set_name + string(spacer_len, ' ') + operation_str;
			summary_string += summarise_intervals(opt, transformed_interval_set_name, max_interval_set_name_len,
													intervals_per_contig, genome_contig_sizes, 1);
			if (opt.operation == eSUMMARISE)
				continue;

			// 
			// write to output file
			// 
			if (opt.output_file_name.length())
			{
				string output_file_name = get_file_name_substitutions(	opt.output_file_name,
																		input_file_name, custom);
				path output_path = path(output_file_name).branch_path();
				if (!output_path.parent_path().empty())
				{
					if (!exists(output_path))
						create_directories(output_path);
					else if (!is_directory(output_path))
						throw std::logic_error("Can't open " + output_file_name + ". " + 
												output_path.string() + " is not a directory.");
	
				}
				ofstream f_output(output_file_name.c_str());
				if (!f_output)
					throw std::logic_error("Can't open " + output_file_name);
				opt.both(3) << "Output= " << output_file_name << "\n\n\n";
				write_intervals(f_output, intervals_per_contig, opt.zero_based_output);
			}
		}

		summary_string += string (80, '*') + "\n\n";
		// have all summaries together once more
		if (input_file_names.size() > 1)
		{
			opt.log() << summary_string;
			if (opt.get_verbosity() > 1)
				opt.verbose() << summary_string;
		}
		opt.summary() << summary_string;


		//  Finished!
		// 
		opt.both() << "Finished!\n\n";

	}
	catch(const std::exception& e)
	{
		// indent error string
		string error_string = e.what();
		replace_all(error_string, "\n", "\n\t");
		std::cerr
				<< "\n"<< string(60, '=') << "\n"
				<< "\nError:\n\t" << error_string 
				<< "\n\n" << string(60, '=') << "\n\n";
		return -1;
	}
	return 0;
}


